#ifndef LINEAGEMCMC_H
#define LINEAGEMCMC_H

#include <math.h>
#include <vector>
#include <iomanip>
#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <stdlib.h>
#include <stdio.h>
#include <sstream>
#include "tree.h"
#include "lineage.h"

#include <gsl/gsl_randist.h>

gsl_rng * r3;

using namespace::std;

// this structure keeps the loglikelihood values for efficient updates
// overload output operator below

 
struct logLikelihood
{
	double logPoolPrior, logLambdaPrior, logTreePrior, logDataLikelihood, logConeLikelihood, logLineageCond, logStateCond, logXiPrior, logEtaPrior;

	logLikelihood()
	{
		logPoolPrior = -1e20, logLambdaPrior = -1e20, logTreePrior = -1e20, logDataLikelihood = -1e20, logStateCond = -1e20, logLineageCond = -1e20, logXiPrior = -1e20; logEtaPrior = -1e20;
	}
	void operator=(logLikelihood rhs)
	{
		logPoolPrior = rhs.logPoolPrior;
		logTreePrior = rhs.logTreePrior; 	

		logDataLikelihood = rhs.logDataLikelihood;

		logLineageCond = rhs.logLineageCond;
		logStateCond   = rhs.logStateCond;		

		logXiPrior = rhs.logXiPrior;
		logLambdaPrior = rhs.logLambdaPrior;
		logEtaPrior = rhs.logEtaPrior;
	}
};

// overloaded output for logLikelihood struct
ostream& operator <<(ostream& os, const logLikelihood s)
{  
	cout<<"\t Log of pool prior: "<<s.logPoolPrior<<endl
	    <<"\t Log of tree prior: "<<s.logTreePrior<<endl
	    <<"\t Log of lambda prior:"<<s.logLambdaPrior<<endl
	    <<"\t Log of xi prior: "<<s.logXiPrior<<endl
	    <<"\t Log of eta prior: "<<s.logEtaPrior<<endl
	    <<"\t Log of conditional lineage probability: "<<s.logLineageCond<<endl
	    <<"\t Log of conditional state probability: "<<s.logStateCond<<endl
	    <<"\t Log of data likelihood: "<<s.logDataLikelihood<<endl;
}


class lineageMCMC
{
	public:
//		basic objects - constructors, declarations, et c.
		lineageMCMC() {r3 = gsl_rng_alloc(gsl_rng_rand48); long seed; seed = rand()*time (NULL) * getpid();  gsl_rng_set (r3, seed);};
		~lineageMCMC() {};
		lineage currentLineage;
		lineage proposedLineage;
	
		vector < vector < vector <int> > > dataD; // data stuctre

		logLikelihood current;
		logLikelihood proposed;

		// Initialization
		
		// MCMC chain updates
		double updateS(int);
		double updateT();
		double updateR(double);

		double updateXi();
		double updateLambda(); // update lambda
		double updateEta();    // update eta
	
		int MetropolisHastings(int, int, string);
	
		// Likelihood, probabilities
		void calculateLogLikelihood(lineage *, logLikelihood *);
		double calculateLogPoolPriorProb(lineage *);
		double calculateLogTreePriorProb(lineage *);

		double calculateLogLineageCondProb(lineage *);
		double calculateLogStateCondProb(lineage *);

		double calculateLogDataLikelihood(lineage *);

		double calculateLogXiPrior(lineage *);
		double calculateLogLambdaPrior(lineage *);				
		double calculateLogEtaPrior(lineage *);

		double recursionProb(lineage*, int, int, int);
		double getPosterior();
		void initializeModel(string, int);
	private:
		int isSiteVariable(int,lineage *);
		double logFactorial(int);
		int vectorRandomElement(vector<double> probVec);
		void normalizeVector(vector < double > & vec);
		int isAllEqual(lineage *, double);
		char drawFunction(int);
		string ToString(int);
		double max(double,double);
		vector <int> whichDiff(lineage * , lineage * );
		double cdf_exp(double,double);
		double matrixExp(double xi, double t, int start, int end);
		double glog(double);
};


double lineageMCMC::max(double a, double b)
{
	if (a >= b)
	{	return(a);	}
	else
	{	return(b);	}
}

string lineageMCMC::ToString(int value)
{
	stringstream ss;
	ss << value;
	return ss.str();
}

double lineageMCMC::calculateLogEtaPrior(lineage * lin)
{
	return(gsl_sf_log(gsl_ran_beta_pdf(lin->eta,1,1)));
}

double lineageMCMC::calculateLogXiPrior(lineage * lin)
{
	double xi = lin-> xi, value;

	value = gsl_ran_exponential_pdf(xi,1);
	
	return(glog(value));
}

double lineageMCMC::calculateLogLambdaPrior(lineage * lin)
{
	return(glog(gsl_ran_beta_pdf(lin->lambda,1,1)));
}

double lineageMCMC::updateLambda()
{
	double logQR = 0, logPR = 0, logLR = 0;

	proposedLineage = currentLineage;
	proposed = current;

	proposedLineage.proposeLambda();

	proposed.logLambdaPrior = calculateLogLambdaPrior(&proposedLineage);
	proposed.logStateCond = calculateLogStateCondProb(&proposedLineage);
//	proposed.logLineageCond = calculateLogLineageCondProb(&proposedLineage);

	logQR = proposedLineage.logLineageProposal;
	logPR = proposed.logLambdaPrior - current.logLambdaPrior;
	logLR = proposed.logStateCond - current.logStateCond;
//	logLR += proposed.logLineageCond - current.logLineageCond;

//	cout<<logPR<<" "<<logLR<<" "<<logQR<<endl;
	return( logPR + logLR + logQR);
}

double lineageMCMC::updateEta()
{
	double logQR = 0, logPR = 0, logLR = 0;

	proposedLineage = currentLineage;
	proposed = current;

	proposedLineage.proposeEta();

	proposed.logEtaPrior = calculateLogEtaPrior(&proposedLineage);
	proposed.logDataLikelihood = calculateLogDataLikelihood(&proposedLineage);
//	proposed.logStateCond = calculateLogStateCondProb(&proposedLineage);

	logQR = proposedLineage.logLineageProposal;
	logPR = proposed.logEtaPrior - current.logEtaPrior;
	logLR = proposed.logDataLikelihood - current.logDataLikelihood;

	return(logQR+logLR+logPR);
}

double lineageMCMC::updateXi()
{
	double logLR = 0, logPR = 0, logQR =0;
	proposedLineage = currentLineage;
	proposed = current;

	proposedLineage.proposeXi();
	proposed.logLineageCond = calculateLogLineageCondProb(&proposedLineage);
	proposed.logStateCond = calculateLogStateCondProb(&proposedLineage);
	proposed.logXiPrior = calculateLogXiPrior(&proposedLineage);
	proposed.logTreePrior = calculateLogTreePriorProb(&proposedLineage);

	logQR = proposedLineage.logLineageProposal;
	logPR = proposed.logXiPrior - current.logXiPrior;
	logPR += proposed.logTreePrior - current.logTreePrior;
	logLR = proposed.logLineageCond - current.logLineageCond;

	return( logPR + logLR + logQR);
}

double lineageMCMC::updateS(int iter)
{
	int i = floor(gsl_ran_flat(r3,0,currentLineage.numPools-1e-10));
	double proposedValue = 0, currentValue = 0, logLR = 0, logQR = 0;

	proposedLineage = currentLineage;
	proposed = current;

//	if (iter >0)
//	{
	proposedLineage.proposePool(i);//dataD,i,iter);
	current.logDataLikelihood = calculateLogDataLikelihood(&currentLineage);
	proposed.logDataLikelihood = calculateLogDataLikelihood(&proposedLineage);
	logLR = proposed.logDataLikelihood - current.logDataLikelihood;

	logQR = 0;//proposedLineage.logLineageProposal;
	
	return logLR+logQR;
}


double lineageMCMC::updateT()
{
	double logQR = 0, logPR =0, logLR =0, uni;

	vector <int> diff;

	uni = gsl_ran_flat(r2,0,1);

	proposedLineage.proposeTree();

	proposed.logLineageCond = calculateLogLineageCondProb(&proposedLineage);
	proposed.logStateCond = calculateLogStateCondProb(&proposedLineage);
	proposed.logDataLikelihood = calculateLogDataLikelihood(&proposedLineage);

	proposed.logTreePrior = calculateLogTreePriorProb(&proposedLineage);
	current.logTreePrior = calculateLogTreePriorProb(&currentLineage);		

	logLR = proposed.logLineageCond - current.logLineageCond;
	logLR += proposed.logStateCond - current.logStateCond;
	logPR += proposed.logTreePrior - current.logTreePrior; 
	logQR = 0;//proposedLineage.roposal;

//	cout<<"T:\t"<<logLR<<"\t"<<logPR<<endl;
	return logQR + logPR + logLR;		
}


double lineageMCMC::updateR(double some)
{
	double logLR, logQR, logPR=0, total;

	int k, K = currentLineage.numLineage,j,M=currentLineage.numSNPs;

	proposedLineage = currentLineage;
	proposed = current;

	//cout<<"Here?"<<endl;
	
	proposedLineage.proposeLineage(some, dataD);
	
	proposed.logDataLikelihood = calculateLogDataLikelihood(&proposedLineage);
	proposed.logLineageCond    = calculateLogLineageCondProb(&proposedLineage);
	proposed.logStateCond      = calculateLogStateCondProb(&proposedLineage);
//	proposed.logLambdaPrior    = calculateLogLambdaPrior(&proposedLineage);

	logLR = proposed.logDataLikelihood - current.logDataLikelihood;
	logQR = proposedLineage.logLineageProposal;
	logPR = proposed.logLineageCond - current.logLineageCond;
	logPR += proposed.logStateCond - current.logStateCond;
//	logPR += proposed.logLambdaPrior - current.logLambdaPrior;

//	cout<<logLR<<"\t"<<logPR<<"\t"<<logQR<<endl;
	return logLR + logPR + logQR;
}


void lineageMCMC::initializeModel(string file, int numLineage)
{
	int j=0;

	currentLineage.setNumLineages(numLineage);
	currentLineage.readData(file, dataD);
	cout<<"The data file has been read."<<endl;
	currentLineage.initializeParameters(dataD);
	cout<<"Parameters randomly initialized."<<endl;
//	cout<<"Lineage model stats:"<<endl<<currentLineage<<endl;
	cout<<"Partitions likelihod calculated."<<endl;
	calculateLogLikelihood(&currentLineage, &current);
	cout<<"Likelihood calculated: "<<current<<endl<<endl;
}


void lineageMCMC::calculateLogLikelihood(lineage * lin, logLikelihood * like)
{
	like -> logPoolPrior = calculateLogPoolPriorProb(lin);
	like -> logTreePrior = calculateLogTreePriorProb(lin);

	like -> logXiPrior = calculateLogXiPrior(lin);
	like -> logEtaPrior = calculateLogEtaPrior(lin);
	like -> logLambdaPrior = calculateLogLambdaPrior(lin);

	like -> logLineageCond  = calculateLogLineageCondProb(lin);
	like -> logStateCond    = calculateLogStateCondProb(lin);

	like ->logDataLikelihood = calculateLogDataLikelihood(lin);
}


double lineageMCMC::calculateLogStateCondProb(lineage * lin)
{
	int j, M = lin->numSNPs, count = 0;

	double value;
	for (j = 0; j < M; j++)
	{
		count += (1-isSiteVariable(j,lin));	
	}

//	cout<<count<<endl;
	value = glog(gsl_ran_binomial_pdf(count, lin->lambda, M));
//	cout<<value<<endl;
	value += double(count)*glog(0.5);

//	cout<<value<<" "<<count<<" "<<lin->lambda<<endl;
	return(value);
}


double lineageMCMC::calculateLogLineageCondProb(lineage * lin)
{
	// felsenstein likelihood
	int M = lin->numSNPs, j = 0, k, K = lin->numLineage;

	double logProb = 0,value=0, value_r, value_l, Z;
	
	vector < double > store;

	store.clear();
	store.resize(K,0);

	for (k = 0; k < K; k++)
	{	
		store[k] = lin->lineageR[k][0];
		lin->lineageR[k][0] = 0;
	}

	Z  = 0.5*recursionProb(lin,0,0,0);
	Z += 0.5*recursionProb(lin,0,0,1);

	Z = Z*2;

	for (k = 0; k < K; k++)
	{	lin->lineageR[k][0] = store[k];	}

	for (j = 0; j < M; j++)
	{
		if (isSiteVariable(j,lin)==1)
		{
			value = 0.5*recursionProb(lin, 0, j, 0);
			value += 0.5*recursionProb(lin, 0, j, 1);
			
			value = value/(1-Z);

			logProb += glog(value);
		}
	}

	return(logProb);
}

double lineageMCMC::cdf_exp(double x, double mu)
{
	return(1-exp(-1*mu*x));
}


double lineageMCMC::recursionProb(lineage *lin, int node, int j, int state)
{
	int k = 0, r_node, l_node, i;
	double xi, r_branch, l_branch,value = 0, value_l=0, value_r=0;

	if (lin->treeT.isTip(node))
	{
		k = lin->treeT.lineageNodeHash[node];
		if (state == lin->lineageR[k][j])
		{	value = 1;	}
		else	
		{	value = 0;	}
	}
	else
	{
		xi = lin->xi;

		r_node = lin->treeT.treeNodes[node].right;
		l_node = lin->treeT.treeNodes[node].left;

		r_branch = lin->treeT.treeNodes[node].rightBranch;
		l_branch = lin->treeT.treeNodes[node].leftBranch;

		value_l = 0;
		value_r = 0;

		for (i = 0; i < 2; i++)
		{

			value_r += matrixExp(xi,r_branch,state,i)*recursionProb(lin,r_node,j,i);
			value_l += matrixExp(xi,l_branch,state,i)*recursionProb(lin,l_node,j,i);

		}
		value = value_r*value_l;
	}
	return(value);
}

double lineageMCMC::matrixExp(double xi, double t, int start, int end)
{
	double value = 0;
	if (start == end)
	{	value = 0.5 + 0.5*exp(-2*xi*t);	}
	else
	{	value = 0.5 - 0.5*exp(-2*xi*t);	}
	return(value);
}

int lineageMCMC::isSiteVariable(int j,lineage * lin)
{
	int K = lin->numLineage, k, total = 0;

	for (k = 0; k < K; k++)
	{
		total += lin->lineageR[k][j];
	}
	if ((total==K)|(total==0))
	{
		return(0);
	}
	else
	{
		return(1);
	}
}

double lineageMCMC::calculateLogDataLikelihood(lineage * lin)
{
	int i, j, k, N = lin->numPools, M = lin->numSNPs, K = lin->numLineage;
	long up =0, down=0, total_n, total_r, dummy;

	double total = 0, p_ij =0;
	double eta = lin->eta;
	total = 0;		
	for (j = 0; j<M; j++)
	{
		for (i = 0;i<N; i++)
		{

			p_ij= 0;
			for (k=0; k<K;k++)
			{	
				p_ij = p_ij + lin->poolS[i][k]*lin->lineageR[k][j];
			}

			p_ij = p_ij - 2*lin->eta*p_ij+eta;

			up = dataD[j][i][1]+dataD[j][i][0];
			down = dataD[j][i][0];

//			cout<<p_ij<<endl;

			if ((p_ij<1-1e-323)&&(p_ij>1e-323))
			{
				total += lin->logBinomial((unsigned)up,(unsigned)down) + (dataD[j][i][0])*gsl_sf_log(p_ij)+(dataD[j][i][1])*gsl_sf_log(1-p_ij);		
			}
			else
			{
				total += -743.75;
			}		
		}
	}

	return total;
}

double lineageMCMC::calculateLogPoolPriorProb(lineage * lin)
{
	int N = lin->numPools, K = lin->numLineage;

	if (isAllEqual(lin, 0.01))
	{  return -N*lin->logFactorial(K-1); }
	else
	{   cout<<"No separate!  No equal!"<<endl;	exit(-1);	}
}


double lineageMCMC::calculateLogTreePriorProb(lineage * lin)
{
	int K = lin->numLineage,  numEvents, i, counter = 0,M = 1;//lin->numSNPs;

	double total = 0, numTips;
	
	vector < int > coalescentNodes;
	vector < double > coalescenceTimes, coalescenceAscend, coalescenceIntervals, coalescentNums;

	total = (K-1)*gsl_sf_log(2)+lin->logFactorial(K-3)-lin->logFactorial(2*K-5) ;

	numEvents = (lin -> treeT.treeNodes.size() - 1)/2;
	numTips = (double)(numEvents + 1);

	coalescenceTimes.clear();  coalescenceTimes.resize(numEvents);
	coalescenceIntervals.clear();  coalescenceIntervals.resize(numEvents);
	coalescentNodes.clear(); 
	coalescentNums.clear(); coalescentNums.resize(numEvents);

	// find all the coalescent times

	for (i = 0; i < lin->treeT.treeNodes.size(); i++)
	{
		if (lin ->treeT.isTip(i) == 0)
		{
			coalescenceTimes[counter] = 1- lin->treeT.findPathAbove(i); 
			coalescentNodes.push_back(i);
			counter ++;		
		}
	}
	coalescenceAscend = coalescenceTimes;

	sort(coalescenceAscend.begin(),coalescenceAscend.end());

	for (i = 0; i < coalescenceTimes.size(); i++)
	{
		if (i == 0)
		{
			coalescenceIntervals[i] = coalescenceAscend[i];
			coalescentNums[i] = numTips - i;
		}
		else
		{
			coalescenceIntervals[i] = coalescenceAscend[i] - coalescenceAscend[i-1];
			coalescentNums[i] = (numTips - i)*(numTips-i-1);
		}
	}

	for (i = 0; i < coalescenceTimes.size()-1;i++)
	{
		total += (gsl_sf_log(coalescentNums[i]/2.0) - lin->xi*coalescenceIntervals[i]*coalescentNums[i]/2.0);
	} 
	
	//cout<<total<<endl;
	return total;
}

int lineageMCMC::isAllEqual(lineage * lin, double tol)
{
	int i, k, check = 1, N = lin->numPools, K = lin->numLineage;
	double total;

	for (i =0; i <N;i++)
	{
		total = 0;
		for (k = 0; k<K;k++)
		{
			total += lin->poolS[i][k];
		}
		if (fabs(total-1)>tol)
		{	check = 0;	}
	}
	return check;
}


double lineageMCMC::logFactorial(int n)
{
	if ((n==1)||(n==0))
	{ 	return 0;	}
	else
	{  	return gsl_sf_lnfact(n);	}
}


void lineageMCMC::normalizeVector(vector<double> & vec)
{
	double total = 0;
	int i =0, check;

	for (i =0; i<vec.size(); i++)
	{	total += vec[i];	}
	
	if (total > 0)
	{
		for (i =0; i<(int)vec.size();i++)
		{
			vec[i] = vec[i]/total;
		}
	}
	else
	{
		for (i =0; i<(int)vec.size();i++)
		{
			vec[i] = 1/(double)vec.size();
		}		
	}
	check = 0;
	for (i = 0; i<vec.size();i++)
	{
		if (vec[i]==0)
		{	check = 1; vec[i] = 1e-20;	}
	}
	if (check == 1)
	{	normalizeVector(vec);	}
	//return vec;
}

int lineageMCMC::vectorRandomElement(vector<double> probVec)
{
	int i = 0, check = 1, counter = 0;
	double  uni;
	vector<double> sumVec;

	normalizeVector(probVec);

	sumVec.clear();	
	sumVec.resize(probVec.size(),0);


	for (i = 0; i<probVec.size(); i++)
	{
		if (i > 0)
		{	sumVec[i] = probVec[i] + sumVec[i-1];}
		else
		{	sumVec[i] = probVec[i];	}
	}
	


	uni = gsl_ran_flat(r2,0,1 - 1e-5);

	while(check)
	{
		if (uni < sumVec[counter])
		{	check = 0;	}		
		else
		{	
			if (counter== probVec.size()-1)
			{	check = 0;	}
			else
			{	counter ++;	}
		}
	}

	return counter;
}

double lineageMCMC::getPosterior()
{
	return(current.logPoolPrior+current.logTreePrior+current.logXiPrior + current.logEtaPrior+current.logLambdaPrior + current.logLineageCond+current.logStateCond + current.logDataLikelihood);
}

char lineageMCMC::drawFunction(int iter)
{
	int call;
	double uni  = gsl_ran_flat(r2,0,1);
	char return_char;
	proposed.logDataLikelihood = calculateLogDataLikelihood(&proposedLineage);
	vector < double > probs;
	probs.clear(); 


	if (iter > -1)
	{
		probs.push_back(4.0); //p 4
		probs.push_back(2.0); //t 2
		probs.push_back(4.0); //r 4
		probs.push_back(0.5); //x

		probs.push_back(1.0); //e 1 
		probs.push_back(1.0); //l 1
	}
	else
	{
		probs.push_back(10.0); //p
		probs.push_back(2.0); //t
		probs.push_back(2.0); //r
		probs.push_back(0.0); //x

		probs.push_back(0.0); //e
		probs.push_back(0.0); //l	
	}

	call = vectorRandomElement(probs);
	
	switch(call)
	{
		case 0:
			return_char = 'p';
		break;
		case 1:
			return_char = 't';
		break;
		case 2:
			return_char = 'r';
		break;
		case 3:
			return_char = 'x';
		break;
		case 4:
			return_char = 'e';
		break;		
		case 5:
			return_char = 'l';
		break;				
	}	
	
//	0.8, 0.8, 
	return return_char;
}

int lineageMCMC::MetropolisHastings(int numIter, int increment_record, string filename)
{
	int i =0,pool, j, k, iter, h, nh;
	char draw;
	double total = 0, uni = 0, alpha;

	string treeOutput;
	vector<int> acceptanceCounts;
	ostringstream oss;
	acceptanceCounts.clear();
	acceptanceCounts.resize(6,0);	

	ofstream fout;
	fout.open(filename.c_str());

	current.logDataLikelihood =  calculateLogDataLikelihood(&currentLineage);		
	current.logLineageCond =  calculateLogLineageCondProb(&currentLineage);
	
	proposed = current;
	proposedLineage = currentLineage;

	cout<<setw(6)<<"Iteration"<<"\t"<<setw(6)<<"logXiPrior"<<"\t"<<setw(6)<<"logEtaPrior"<<"\t"<<setw(8)<<"logTreePrior"<<"\t"<<setw(6)<<"logPoolPrior"<<"\t"<<setw(6)<<"logLLK"<<"\t"<<setw(8)<<"logCondProb"<<"\t"<<setw(6)<<"logStateProb"<<"\t"<<setw(6)<<"logLambda"<<endl;
			
	for (iter = 0; iter < numIter; iter++)
	{	 
		proposedLineage = currentLineage;
		proposed = current;

		draw = drawFunction(iter);
	//	cout<<draw<<endl;
		switch(draw)
		{
			case 't':
				alpha = updateT();
			break;
			case 'x':
				alpha = updateXi();
			break;
			case 'r':
				alpha = updateR(0.05);
			break;
			case 'p':
				alpha = updateS(iter);
			break;
			case 'e':
				alpha = updateEta();
			break;
			case 'l':
				alpha = updateLambda();
			break;

			default:
				cerr<<"No draw"<<endl;
				exit(-1);
			break;
		}	

		alpha = min(alpha, 0.0);		

		uni = gsl_sf_log(gsl_ran_flat(r2,epsilon,1-epsilon));

		if (uni < alpha)
		{	

			current = proposed;
			currentLineage = proposedLineage;

			switch(draw)
			{
				case 't':
					acceptanceCounts[0]++;
				break;
				case 'p':
					acceptanceCounts[1]++;
				break;
				case 'r':
					acceptanceCounts[2]++;
				break;
				case 'x':
					acceptanceCounts[3]++;
				break;
				case 'l':
					acceptanceCounts[4]++;
				break;
				case 'e':
					acceptanceCounts[5]++;
				break;

			}			
		}
		
		if (iter % increment_record == 0)
		{	
			cout<<setw(6)<<iter<<"\t"<<setw(6)<<current.logXiPrior<<"\t"<<setw(6)<<current.logEtaPrior<<"\t"<<setw(8)<<current.logTreePrior<<"\t"<<setw(6)<<current.logPoolPrior<<"\t"<<setw(6)<<current.logDataLikelihood<<"\t"<<setw(8)<<current.logLineageCond<<"\t"<<setw(6)<<current.logStateCond<<"\t"<<setw(6)<<current.logLambdaPrior;
			cout<<setw(6)<<endl;
			
			fout<<"Iter: "<<iter<<endl;
			fout<<"Probs: "<<current.logXiPrior<<" "current.logEtaPrior<<" "<<current.logTreePrior<<" "<<current.logPoolPrior<<" "<<current.logDataLikelihood<<" "<<current.logLineageCond<<" "<<current.logStateCond<<" "<<current.logLambdaPrior<<" "<<endl;
			fout<<"Number of lineages: "<<currentLineage.numLineage<<endl;
			fout<<'#'<<endl<<"Pools:"<<endl;
			for (i = 0; i<currentLineage.numPools; i++)
			{
				for (k = 0; k< currentLineage.numLineage; k++)
				{
					fout<<currentLineage.poolS[i][k]<<" ";
				}
				fout<<endl;
			}
			fout<<endl<<'#'<<endl<<"Lineage:"<<endl;		
			for (k = 0; k<currentLineage.numLineage; k++)
			{
				for (j = 0; j< currentLineage.numSNPs; j++)
				{
					fout<<currentLineage.lineageR[k][j];
					if (j < currentLineage.numSNPs - 1)
					{	fout<<" ";	}
				}
				fout<<"\n";
			}
			fout<<'#'<<endl<<"Xi:"<<endl;
			fout<<currentLineage.xi<<endl;
			fout<<'#'<<endl<<"Eta:"<<endl;
			fout<<currentLineage.eta<<endl;
			fout<<'#'<<endl<<"Lambda:"<<endl;
			fout<<currentLineage.lambda<<endl;
			fout<<'#'<<endl<<"Newick:"<<endl;
			fout<<currentLineage.treeT.outputTree()<<endl;
			fout<<'#'<<endl<<"Tree:"<<endl;
			fout<<currentLineage.createLineageGraph()<<endl;
			fout<<'#'<<endl<<"Special"<<endl;
			fout<<currentLineage.treeT<<endl;		
			fout<<'#'<<endl;
		}
	}

	cout<<endl<<endl<<"Acceptance counts:";
	for (i = 0; i<acceptanceCounts.size(); i++)
	{
		cout<<acceptanceCounts[i]<<" ";
	}
	cout<<currentLineage<<endl;
	currentLineage.treeT.displayTree();
	return 1;
}

double lineageMCMC::glog(double value)
{
	if (value > 1e-323)
	{
		return(gsl_sf_log(value));
	}
	else
	{
		return(-743.75);
	}
}

#endif /* LINEAGEMCMC_H */
